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We discuss recently introduced numerical linked-cluster (NLC) algorithms that allow one to ob- 
tain temperature-dependent properties of quantum lattice models, in the thermodynamic limit, from 
exact diagonalization of finite clusters. We present studies of thermodynamic observables for spin 
models on square, triangular, and kagome lattices. Results for several choices of clusters and ex- 
trapolations methods, that accelerate the convergence of NLC, are presented. We also include a 
comparison of NLC results with those obtained from exact analytical expressions (where available) , 
high-temperature expansions (HTE), exact diagonalization (ED) of finite periodic systems, and 
quantum Monte Carlo simulations. For many models and properties NLC results are substantially 
more accurate than HTE and ED. 
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I. INTRODUCTION 

Accurate quantitative calculations of finite- 
temperature properties of quantum lattice models 
are a challenging task [l|, One of the few general 
methods that works directly in the thermodynamic 
limit, is that of high-temperature expansions (HTEs), 
where properties of the system are expanded in powers 
of inverse temperature, (3 — 1/T 3]. These expansions, 
carried out to order (3 N (where N is typically around 
10), provide highly accurate numerical results at high 
temperatures. However, below some temperature T s 
related to the relevant microscopic energy scale of 
the system, the series diverges. Such a divergence 
need not be related to any finite-temperature phase 
transitions or long-range correlations. For example, in 
low dimensional or frustrated spin systems, often there 
is either no finite-temperature phase transition or such 
a transition occurs well below the exchange constant J. 
The microscopic energy scale J still controls the radius 
of convergence of the high-temperature series. 

Beyond the radius of convergence [3 > (3 S , series ex- 
trapolation methods [4] allow one to calculate thermody- 
namic properties, but their reliability remains uncertain. 
Our motivation for developing the numerical linked clus- 
ter (NLC) method is to be able to obtain the properties 
of these systems, in the thermodynamic limit, for (3 > (3 S 
in a more reliable way, especially if the correlations in the 
system remain short ranged. NLC uses the linked cluster 
basis of HTE, but replaces the expansion in (3 by an exact 
numerical calculation of the properties of linked clusters 
at any temperature [H, Q . In any practical implementa- 
tion of NLC, only the contributions from clusters up to 
some maximum size are included. This can lead to highly 
accurate properties of the thermodynamic system, even 
beyond the radius of convergence of HTE, provided the 
correlations are short ranged. Thus, NLC helps to sepa- 



rate cases where the failure of HTE is due to its analytic 
structure in the complex plane, from where the correla- 
tions truly exceed the largest clusters studied. We would 
like to emphasize that this does not imply that compara- 
ble results for thermodynamic systems can be obtained 
simply by exact diagonalization (ED) of individual pe- 
riodic clusters of size comparable to the maximum size 
used in NLC. We will show that NLC can be substantially 
more accurate than ED for finite-temperature properties 
of a thermodynamic system. Furthermore, one can ac- 
celerate the convergence of NLC, even when correlation 
length exceeds the largest cluster studied, by using se- 
quence extrapolation techniques, which are in many ways 
analogous to series extrapolation methods [J, 0] ■ 

We present in this paper a detailed exposition of the 
NLC algorithm. The basic method was already presented 
in Ref. [g| . Here we discuss in detail the different choices 
of clusters that can be used to build the numerical expan- 
sion. We also detail different extrapolation methods that, 
like Pade approximants for HTE, allow one to accelerate 
the convergence of NLC. These methods are especially 
relevant for the application of NLC to models in which 
correlations build up rapidly with reduction in temper- 
ature. Comparisons between results obtained by means 
of NLC with known techniques such as exact diagonal- 
ization (ED), quantum Monte Carlo (QMC) simulations, 
and HTE, are also presented. 

The exposition is organized as follows. In Sec. UH we 
introduce the basis of NLC. We then present (Sec. IIII[) 
an overview of different sequence extrapolation methods 
that can help accelerate the convergence of NLC. The 
rest of the manuscript is devoted to show how to build 
the series for spin models on square (Sec. IIV[) . triangular 
(Sec. [VJ, and kagome (Sec.lVI[) lattices. Different choices 
of clusters are discussed in detail, and applied to Ising, 
XY, and Heisenberg models. Finally, the conclusions are 
presented in Sec. IVIH 
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II. BASIS OF NLC 



The fundamental basis for a linked cluster expansion, 
for some extensive property P of an infinite lattice C, is 
the relation [j| Q 



P{L)/N = Y,L{c) xWp(c), 



(1) 



where the left hand side is the value of the property P 
per lattice site in the thermodynamic limit. On the right- 
hand side L(c) is the so-called lattice constant, which 
is the number of embeddings of the cluster c, per lat- 
tice site, in the lattice C (explicit examples will be given 
later). Wp(c) is the weight of the cluster c for the prop- 
erty P. The latter is defined recursively by the principle 
of inclusion and exclusion Q , 



W P (c) = P(c) -J2Wp(s), 



(2) 



where P(c) is the property P calculated for the finite 
cluster c. The sum on s runs over all subclustcrs of c. In 
HTE, for every cluster, P(c) and hence its weight Wp(c) 
are expanded in powers of (3 and only a finite number of 
terms are retained. In NLC an exact diagonalization of 
the cluster is used to calculate P(c) and Wp(c) at any 
temperature. [Notice that in contrast to the clusters used 
in ED studies, the ones in Eq. ([2]) do not have periodic 
boundary conditions.] The exact numerical calculation 
of P(c) allows NLC to build more bare information of 
the system than HTE. 

There is another aspect in which the NLC scheme is 
fundamentally different from HTE, and that can be used 
to ones advantage. In HTE, the choice of clusters is dic- 
tated solely by the need to get the power series expansion 
in [3 to as high an order as possible. This typically means 
that clusters are defined and ordered by the number of 
bonds. In NLC, one has substantial freedom to select 
the set of clusters and the order in which they are con- 
sidered. One can arrange them by number of sites, by 
number of bonds or, as we will see below, one can even 
consider a reduced set of clusters. In order to obtain 
correct results at high temperatures, one requirement is 
that with increasing order, the cluster weights, when ex- 
panded in inverse temperature, should give the correct 
HTE coefficients as well. This ensures that when HTE 
converges, NLC gives results that are identical to HTE. 
However, as we will see, NLC expansions may involve a 
choice of clusters that sacrifice efficiency in the order to 
which they give the exact HTE coefficients, for better 
results at intermediate and low temperatures. 

As discussed before, when correlations are shorter 
ranged than the size of clusters that are included in 
any implementation of NLC, the NLC results are accu- 
rate without any need for extrapolations. On the other 
hand, for systems with an ordered ground state, such 
an implementation must eventually break down as the 
temperature is lowered (larger clusters will start making 



large contributions). In the next section we discuss some 
"tricks" that can be used to accelerate the convergence 
of the direct sum in Eq. (p} . This can lead to accurate 
thermodynamic results at temperatures where correla- 
tion length far exceeds the sizes of the cluster included 
in the sum. 



III. SEQUENCE EXTRAPOLATIONS 

In this section, we consider the general problem of es- 
timating P(C) in Eq. {]}, when the weights Wp are only 
known for clusters up to a given size. In order to produce 
a sensible extrapolation scheme one can group clusters 
together and define 



S n =J2 L (cn) x W P (c n ), 



(3) 



where all clusters c n share a given characteristic, which 
could be, for example, the number of bonds, sites, etc. 
Equation can be rewritten as 



P(C)=J2Sn, 

n 

and one can define partial sums as 

71 

P n (£) = J2Si 



(4) 



(5) 



So, our goal is to estimate P{£) = lim, woo P n from a 
sequence {Pn}, which is known for n = 1, . . . , N. 

Several methods have been developed to accelerate the 
convergence of such sequences. An extensive review with 
references to original papers can be found in Ref . [ij , and 
they are similar to series extrapolation methods. We will 
briefly describe here three methods that we have imple- 
mented, and that have proven to be particularly useful 
in accelerating the convergence of NLC. These methods 
are known as the Wynn's (e) algorithm |4], the Brezin- 
ski's (#) algorithm and the Euler's transformation 
Q . A very important topic that is discussed neither here 
nor m Refs. @, is one of error estimation. This is 
because most studies of errors associated with such ex- 
trapolation methods depend on the underlying function 
satisfying certain properties. As one might expect, these 
properties, in general, cannot be verified for many of the 
problems one finds in statistical mechanics. 

Wynn's algorithm is defined as follows [37j 



,0) 



e (0) - P 



Jk-2) 
-n+1 



Ae 



(fe-i) ■ 



(6) 



where the discrete differentiating operator A is only ap- 
plied to subscripts 



Ae 



(fc-i) _ Jfc-i) _ jk-i) 



(7) 
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Within this scheme the even entries el 2 ^ are expected to 
converge to P{C), while the odd ones e^ 2fc+1 ^ diverge. 
Nonlinear sequence extrapolations usually display this 
behavior, and it implies that one has to be careful with 
round-off errors. 

One should notice that two iterations are needed for 
each level of improvement so the new sequence is two 
terms shorter. We refer to each level of improvement as 
a cycle. For the problems we have studied so far by means 
of NLC, Wynn's algorithm has been, in general, the most 
successful in extending the region of convergency to lower 
temperatures. 

Brezinski's algorithm is defined as follows: 



{ n 1] = 0, 6 
(2fc+l) _ g(2k-l) 



P, 



(2k) A/ n(2fc+l) 



42k) At Wi At Wi 



'n+l 



(8) 



once again the discrete differentiating operator A is only 
applied to subscripts as in Eq. ((TJ and 



(9) 



As for Wynn's algorithm, only even entries are ex- 
pected to converge to P(C). Once again, we have a cycle 
of improvement after two iterations, and for the Brczin- 
ski algorithm three terms are lost in each cycle. This 
fact, together with the second derivative in the denomi- 
nator in #i 2fc+2 - ) [Eq. ([5])] , reduce the number of cycles one 
can perform with the Brezinski algorithm as compared to 
Wynn's algorithm. 

Finally, for alternating series, i.e., series in which terms 
S n [Eq. ((4])] alternate in sign [S n = (— l)"it n ], there is 
a powerful tool known as the Euler transformation Q. 
With this algorithm P^ (£) is approached by the sum (n 
is even) 



u Q - Ul + u 2 



U n -1 



E 



2s+i L '• 



(10) 



where A is the forward difference operator 

Au n = u n+ i - u n , 
A 2 u n = u n+2 - 2u n+ i + u n , 
A 3 u n = u n+3 - 3u„ +2 + 3w„+i - it, 



(11) 



It is always advisable to do the sum of a small number of 
terms directly, through term n — 1, and then apply the 
transformation. 

As we will show later, we have found Euler's trans- 
formation to be particularly useful for calculations of 
the specific heat. Having stated that NLC provides a 
scheme that similar to HTE allows for systematic ex- 
trapolations, we should add some remarks here. What 
the NLC scheme misses is the analytic structure of HTE. 



These have proven particularly useful in studies of criti- 
cal phenomena, where the region of diverging correlation 
length has also been addressed using Pade extrapolations 
[|| . We have not yet investigated if NLC can be used to 
study critical phenomena, as our focus has been on mod- 
els that do not order down to very low temperatures. 

The analytic structure of HTE also allows for changes 
of variables, and in some cases this can be used in very 
ingenious ways. For example, recent work by Bernu and 
Misguich Q , has shown that by converting the expansion 
for entropy in the inverse-temperature variable to an ex- 
pansion for entropy in the internal energy, one can de- 
velop a powerful extrapolation scheme that can build the 
ground-state energy and low-temperature power-law be- 
havior into the extrapolation of high-temperature expan- 
sions. However, we note that such a scheme is only known 
for the entropy (and related quantities) and not for sus- 
ceptibilities or correlation functions. We will show that 
for the specific heat (and entropy) of two-dimensional 
Heisenberg antiferromagnets the NLC scheme works bet- 
ter than direct extrapolation of HTE, but the method of 
Bernu and Misguich is better (at least for square and 
triangular lattice Heisenberg models) in that it allows 
estimates all the way down to T = 0. 

In what follows we study thermodynamic quantities 
(internal energy, entropy, specific heat, and uniform sus- 
ceptibility) of spin models to show how different clusters 
and extrapolation techniques can be used to build NLC 
on square, triangular, and kagome lattices. 



IV. SQUARE LATTICE 

In this section we consider the square lattice. We dis- 
cuss three different cluster schemes that we have used to 
build our NLC expansions. 

The first, and most natural, choice is to consider all 
clusters and order them by the number of bonds. This se- 
lection has been called "Weak Embeddings" in the series 
expansion literature Q, and is typically used for HTE. In 
Fig. [T] we show all clusters that have up to three bonds, 
and their lattice constant. 

Notice that one must include the single site cluster, 
which corresponds to zero bonds. It dominates observ- 
ables such as the entropy at very high temperatures. For 
the smallest bond clusters, such as the one with one bond 
(c = 2) and the first with two bonds (c = 3), it is easy to 
realize that L(c) = 2 since in the square lattice they can 
be only placed horizontally, and vertically. The second 
cluster with two bonds (c = 4) can be placed in four ways 
[L(c) = 4], as one can realize rotating it by 90°, and so 
on. 

It is convenient to group together all clusters that have 
the same Hamiltonian, and diagonalize them just once. 
Since, the Hamiltonian depends on the topology of how 
the sites are connected, we call them topological clusters. 
Looking at our example, clusters c = 3 and 4, or c = 6-9 
have the same topology. For calculating thermodynamic 
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FIG. 1: All clusters with up to three bonds and their lattice 
constant for the square lattice. 



properties, all clusters with the same topology make the 
same contribution. For the square lattice, we have cal- 
culated all possible clusters, and their lattice constants, 
up to 14 bonds. The number of topological clusters, and 
sum of L(c), when grouped by their number of bonds is 
presented in Table. |TJ 

A second choice is to identify clusters by sites. When 
building the Hamiltonian for such expansion, one places 
all possible bonds that connect any pair of sites in the 
graph. This leads to a set of clusters and embeddings 
that are called "strong embeddings" in the series expan- 
sion literature [|[ , and is typically used for generating the 



TABLE I: Number of topological clusters and sum of the lat- 
tice constants for clusters with up to 14 bonds in the square 
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"low-temperature expansions" for Ising models. They 
can be used to generate HTEs as well. However, differ- 
ent clusters, with a given number of sites, will contribute 
to HTE in different orders. Thus, the order to which 
the HTE is correct will be determined by a subset of the 
clusters with the same number of sites that happen to 
contribute in the lowest order. On the other hand, hav- 
ing lots of compact clusters, with multiple connectivity 
between points, could mean that they can give better re- 
sults beyond the radius of convergence of HTE. In Fig. [2] 
we show all clusters that have up to four sites, and their 
lattice constant. 
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FIG. 2: All clusters with up to four sites and their lattice 
constant for the square lattice. 

By comparing Figs. [1] and [5] one can see that the latter 
never includes graphs such as c = 7 of the former one, 
i.e., all squares are always closed in the site expansion, 
hence the name "strong embeddings" . In addition, while 
each site in the square lattice has four nearest-neighbor 
sites, each bond has six nearest-neighbor bonds, which 
implies that the number of bond clusters increases much 
faster than the number of site clusters. In the latter case, 
we have calculated all possible site graphs with up to 16 
sites. Their number of topological clusters and sum of 
lattice constants, when grouped by number of sites, is 
shown in Table [J 

Looking at Tables U and |TT] one can see that for NLC 
calculations of bond and site based expansions the main 
limitation is the computing time (too many clusters) [38j ] , 
and not the memory as is usual for full diagonalization 
studies of clusters with periodic boundary conditions. 
Within NLC one can, however, change that order of lim- 
itations considering more complicated (larger) building 
units for the clusters. Hence, drastically reducing the 
number of different clusters to consider Q. 

A natural selection of a larger building unit in the 
square lattice is, of course, the elementary plaquette or 
square. In this consistent NLC scheme requires 
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TABLE II: Number of topological clusters and sum of lattice 
constants for clusters with up to 16 sites in the square lattice. 

No. of sites No. of topological clusters ^( c ) 



1 11 

2 1 2 

3 1 6 

4 3 19 

5 4 63 

6 10 216 

7 19 760 

8 51 2725 

9 112 9910 

10 300 36 446 

11 746 135 268 

12 2042 505 861 

13 5450 1903 890 

14 15 197 7 204 874 

15 42 192 27 394 666 

16 119 561 104 592 937 



that each bond belongs to only one square. This means 
that we build our clusters out of every alternate square. 
Different squares can only share sites, which are the ze- 
roth order of the square expansion, and are properly sub- 
tracted when calculating the weights in Eq. @. In Fig. 
[3] we show all clusters, with up to three such squares, 
required for a consistent square based NLC expansion. 

The calculation of all possible clusters up to five 
squares (up to 16 sites) is in this case very simple. In 
Table Mil we show the results for the number of topolog- 
ical clusters and sum of their lattice constants. 

In the next subsections we apply the different NLC ex- 
pansions detailed above to several well known spin mod- 
els. All calculations were done on (3.2 GHz) Pentium IV 
personal computers in times that span between 16 h for 
the square based NLC expansion (up to 5 squares) and 

c L(c) 
• 1 1 

I~I 2 1/2 

rf 1 3 1 




FIG. 3: All topological clusters with up to three squares and 
their lattice constant for a square expansion. 



TABLE III: Number of topological clusters and sum of the 
lattice constants for clusters with up to five squares in the 
square lattice. The cluster with zero squares is the single site 
graph. 

No. of squares No. of topological clusters L(c) 
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1 
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1/2 
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1 
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19/2 
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11 


63/2 



60 h for the bond based NLC expansion (up to 13 bonds). 

A. Heisenberg model 

We now consider the antiferromagnetic Heisenberg 
model (AFHM) on the square lattice. Its Hamiltonian 
can be written as 

H = ^S J -S„ (12) 

where we have chosen the coupling constant to be unity, 
and the sum runs over nearest-neighbor spins. 

The AFHM on the square lattice is known to have an 
ordered ground state with long-range antiferromagnetic 
correlations [1(3] • This model can be efficiently simulated 
using QMC techniques, such as stochastic series expan- 
sions [ll|. QMC methods enable one to study much 
larger system sizes than the ones accessible with exact 
diagonalization, although the classes of models that can 
be addressed are limited by the sign problem [H, QJ, El . 

We start by studying the temperature dependence of 
the AFHM energy. In Fig. 01 we show a comparison of the 
bare sums for the bond, site, and square NLC expansions, 
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FIG. 4: (Color online) Energy as a function of temperature for 
the antiferromagnetic Heisenberg model on the square lattice. 
Bare NLC sums are compared with QMC results for a 100 x 
100 lattice. 
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with QMC results using the SSE technique [l5|. Several 
issues are apparent in Fig. 2J (i) All NLC expansions 
give the correct result at high temperatures, (ii) Direct 
NLC sums can converge below T = 1, something that is 
not possible using HTE. (hi) The site expansion, which is 
closer in spirit to low-temperature expansions, performs 
better than the bond expansion (closer to HTEs). This 
occurs even though the site expansion, up to the same or- 
der, is less demanding computationally than the bond ex- 
pansion [there are many topological clusters in the bond 
expansion (1844) with 13 bonds and 14 sites, while in the 
site expansion only clusters with up to 13 sites were diag- 
onalized] . (iv) The direct sum of the largest size cluster 
expansion, the square expansion in this case, converges 
to the lowest temperature (T ~ 0.5), to be compared 
with (T ~ 0.8) for the site expansion. 

The AFHM on the square lattice is known to develop 
antiferromagnetic correlations at relatively high temper- 
atures, i.e., it is the kind of model for which the direct 
NLC sum cannot converge up to very low temperatures. 
Now, the natural question that arises is how low in tem- 
perature can one go by using the sequence extrapola- 
tion techniques discussed in Sec. IIIII In Fig. [5] we show 
two such extrapolations compared with the QMC results. 
The subindex following the name of each extrapolation 
stands for the numbers of cycles of improvements done. 
In addition, for each cycle of improvement there are, in 
general, several terms. In what follows, when not explic- 
itly specified otherwise, we will be showing the highest 
one. 
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FIG. 5: (Color online) Energy as a function of temperature 
for the antiferromagnetic Heisenberg model on the square lat- 
tice. Different extrapolations for the NLC site expansion are 
compared with QMC results for a 100 x 100 lattice, and with 
exact diagonalization results for a 4 x 4 cluster (with periodic 
boundary conditions). 

Figure [5] shows that extrapolations can indeed work 
very well within NLC. We only show results for extrapo- 
lations of the site expansion since for this expansion we 
obtain the best results. Already at the level of the bare 
sum we saw that the site expansion performs better than 
the bond expansion. On the other hand, the square ex- 
pansion, which produces the best results for the direct 



sum, has too few terms to allow for a successful extrap- 
olation scheme to work. [One can also realize that up to 
13 sites (with 8739 different topological clusters), the site 
expansion has explored much more of the lattice than the 
square expansion (with only 21 different topological clus- 
ters).] For the site expansion both Wynn's and Brczin- 
ski's algorithms converge, and agree with QMC, down 
to T ~ 0.15. One should notice, however, that while 
Wynn's results are smooth and on top of the QMC re- 
sults all the way down to T ~ 0.15, some points in the 
Brezinski's extrapolation fall away from that curve. This 
is the kind of numerical problem one can run into after 
several orders of extrapolations. However, with the ex- 
ception of these few points, Brezinski's results are still 
well converged and on top of the QMC data. Also shown 
is the exact diagonalization result for a 4 x 4 cluster with 
periodic boundary conditions. It shows substantial finite- 
size effects already above T > 1. 

We now consider other thermodynamic quantities of 
interest, such as the entropy and the specific heat. For 
the former, we have already shown 5] NLC results for 
Brezinski's and Wynn's extrapolations compared with 
the results of Bernu and Misguich [9[ (obtained by inte- 
grating their specific heat curves). They exhibit a perfect 
agreement down to T ~ 0.3, where S ~ 0.05. Here, we 
will show the results for the specific heat. 
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FIG. 6: (Color online) Specific heat as a function of tempera- 
ture for the antiferromagnetic Heisenberg model on the square 
lattice. Different extrapolations for the NLC site expansion 
are compared with the results of Bernu and Misguich [j| and 
with exact diagonalization results of small clusters (with pe- 
riodic boundary conditions). 

In Fig. [5] we compare NLC results for the specific heat 
(after extrapolation) with the ones obtained by Bernu 
and Misguich [9j. Both approaches basically agree in 
the location of the specific heat peak, although they give 
slightly different peak values. Since NLC results have 
not converged fully below the peak, they may lead to 
the same results as Bernu and Misguich if a few more 
orders are included. Still, for C v , NLC performs much 
better than direct Pade extrapolation of HTEs. QMC 
simulations for the specific heat also suffer from large 
errors (C v has to be obtained as derivative of the energy), 
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and do not allow one to select any of the linked cluster 
results over the other [l5|. 

We also show in Fig.[6]the specific heat results from the 
exact diagonalization of 4 x 4 and 3x3 clusters with pe- 
riodic boundary conditions (PBCs), which helps to give 
an idea of the order of finite size effects in this model. 
They lead to a peak in the specific heat that is neither 
correct in its position nor height. There is another point 
to consider when comparing NLC with ED especially in 
dimensions greater than 1. In the former case, we con- 
sider all possible clusters that make the NLC expansion 
consistent, without any biasing for any type of order. In 
exact diagonalization studies, PBCs can bias the system 
towards or away from certain types of order. For exam- 
ple, in the AFHM the 3 x 3 cluster with PBCs has much 
bigger finite-size effects, because PBCs frustrates antifer- 
romagnetism. This issue may become important when 
the model under consideration has several competing or- 
ders, and the selection of a particular finite size cluster 
may artificially favor one order over the other. NLC, sim- 
ilar to HTE, does not suffer from this problem, and gives 
an unbiased answer for the thermodynamic properties. 



B. XY model in a staggered transverse field 

As discussed in Ref. [5(, NLC is ideal to study models 
that stay short ranged at all temperatures or models in 
which correlations build up slowly. In this section we 
discuss an example of the former case, the XY model 
in a staggered transverse field. Its Hamiltonian can be 
written as 

n = Y J (SfS- + SfS]) + A^(-l) 4 * +i ^ 2 , (13) 

where we have chosen the XY coupling to be unity, the 
sum runs over nearest-neighbor spins, and the last 
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FIG. 7: (Color online) Specific heat as a function of tempera- 
ture for the XY model in a staggered transverse field. A = f 
so that the system is in the insulating regime. Direct sums 
and different extrapolations for the NLC site expansion are 
compared with exact diagonalization results for a 4 x 4 cluster 
(with periodic boundary conditions). 



term describes the staggered field with strength A. 

The XY model in a staggered transverse field can be 
mapped onto a hard-core boson model, at half filling, 
with a staggered site-dependent chemical potential. Such 
a model has been recently studied in one [16|, two [171 ] . 
and three dimensions [l8(. In ID, due to its mapping 
to noninteracting fcrmions, one can realize that there is, 
at zero temperature, a phase transition from a super- 
fluid to an insulating phase for A c = 0, i.e., any finite 
A produces an insulating phase [16|. In two dimensions 
this model has been studied by means of QMC simula- 
tions, in this case the transition between the superfluid 
(also Bose-Einstein condensed phase) and the insulating 
phase occurs, at zero temperature, when A c = 0.992 [l~7l ]. 
Finally, in three dimensions this model has been used to 
rigorously prove the existence of Bose-Einstein condensa- 
tion and Mott insulating phases when tuning the strength 
of the staggered chemical potential [HI . 

In what follows we consider the two-dimensional case. 
In Fig. [7] we show the specific heat as a function of the 
temperature in a case where the system is insulating at 
zero temperature. We have chosen A = 1, close to the 
critical value A c for the superfluid-insulator transition. 
Due to the presence of a gap at zero temperature we 
can be certain that correlations stay short ranged at all 
temperatures. However, they are not negligible, and the 
direct sum of the NLC site expansion exhibits a small os- 
cillation below the peak in the specific heat. A fully con- 
verged result, at all temperatures, can be obtained after 
just one cycle of improvement with Wynn's algorithm or 
using Euler's transformation. Figure [7] also shows that 
the exact diagonalization results for a 4 x 4 cluster with 
periodic boundary conditions still suffer from apparent 
finite size effects. 
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FIG. 8: (Color online) Specific heat as a function of temper- 
ature for the XY model in a staggered transverse field, with 
A = 0.5. Direct sums and different extrapolations for the 
NLC site expansion are compared with exact diagonalization 
results for a 4 x 4 cluster (with periodic boundary conditions). 



Once in the regime of A where the system is superfluid 
at zero temperature, the convergence of the NLC direct 
sum does not reach very low temperatures, but sequence 
extrapolations allow one to reach the region below the 
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peak in the specific heat. This can be seen in Fig. [8] 
As expected, in this case the difference with ED is even 
larger than when A > A c . 

C. Ising model 

To conclude this section on spin models on the square 
lattice, we discuss in this subsection the Ising model 

(id) 

which is an exactly soluble classical model (l9j . 

In two dimensions, the Ising model is known to ex- 
hibit a finite-temperature transition between an ordered 
(gapped) phase, and a disordered high-temperature 
phase. As shown in Fig. [9l the specific heat exhibits 
a divergence at the critical point, which is known to be 
logarithmic [l9[. For NLC, the Ising model reduces to a 
counting problem as the Hamiltonian is already diagonal. 
In Fig. [5] we show direct sums for the site based expan- 
sion up to 15 sites. Surprisingly, Wynn's extrapolations 
allow one to obtain very good estimates of the specific 
heat very close to the critical point. 
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FIG. 9: (Color online) Specific heat as a function of temper- 
ature for the Ising model. Direct sums and different extrap- 
olations for the NLC site expansion are compared with exact 
analytical results. 

We should stress, however, that for finite-temperature 
phase transitions with power-law singularities, HTE is 
probably the most efficient way to go. This is because 
such singularities can be built into the extrapolation. We 
do not know if this is possible within NLC. 

V. TRIANGULAR LATTICE 

In this section we consider the triangular lattice. We 
discuss three different choices of basic clusters that allow 
one to build a consistent NLC expansion. 

The first possibility is, as in the square lattice, to build 
all possible clusters with up to a given number of bonds. 



TABLE IV: Number of topological clusters and sum of the 
lattice constants for clusters with up to 12 bonds in the trian- 
gular lattice. The cluster with bonds is the one site graph. 

No. of bonds No. of topological clusters ~}2 L(c) 



1 1 

1 1 3 

2 1 15 

3 3 91 

4 5 603 

5 12 4215 

6 28 30 535 

7 72 226 905 

8 198 1718 454 

9 590 13 207 569 

10 1817 102 707301 

11 5886 806 366139 

12 19 753 2 086 381866 



Their number of topological clusters, and sum of L(c), 
when grouped by the number of bonds is presented in 
Table 13 

From the number of topological clusters and the sum 
of lattice constants in Table llVl one can see that the num- 
ber of graphs in the triangular lattice grows much faster 
than in the square lattice. This is because in the trian- 
gular lattice each bond has ten nearest-neighbor bonds, 
as opposed to six in the square lattice. (The number of 
nearest-neighbor bonds determines the rate of growth of 
the number of possible clusters.) 

A second natural choice is to build all clusters with 
up to a given number of sites. Once again, when build- 
ing the Hamiltonian for such clusters one needs to place 
the maximum number of bonds possible in them. This 



TABLE V: Number of topological clusters and sum of the lat- 
tice constants for clusters with up to 13 sites in the triangular 
lattice. 

No. of sites No. of topological clusters L{c) 



1 


1 


1 


2 


1 


3 


3 


2 


11 


4 


4 


44 


5 


8 


186 


6 


22 


814 


7 


54 


3652 


8 


156 


16 689 


9 


457 


77 359 


10 


1424 


362 671 


11 


4505 


1716 033 


12 


14 791 


8182 213 


13 


49 138 


39 267 086 
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selection of building blocks for the NLC expansion pro- 
vides a significant reduction in the number of clusters one 
needs to consider as compared to the bond expansion. 
(Each site in the triangular lattice has only six nearest 
neighbor sites.) In addition, having more compact clus- 
ters, this expansion performs better and allows for better 
extrapolations in the intermediate and low-temperature 
regimes. The number of topological clusters and sum 
of lattice constants for the site based NLC expansion is 
shown in Table fVl 

A triangular lattice is made out of triangles, so it is 
also possible to develop a NLC for the triangular lattice 
where the clusters consist of closed triangles. However, 
in this consistent NLC scheme requires that one 

restricts the calculation to a single site plus clusters of 
up (or down) pointing triangles only. The reason for this 
restriction is that all bonds of the triangular lattice be- 
long to a unique up (or down) pointing triangle. Different 
triangles should only share sites. The number of topo- 
logically distinct clusters with 1 through 8 triangles, and 
the sum of their lattice constant is shown in Tabic IVII 
(We have grouped them by the number of triangles.) 

TABLE VI: Triangular lattice number of topological clusters 
and sum of the lattice constants for clusters with up to eight 
triangles. The cluster with triangles is the single site. 

No. of triangles No. of topological clusters X^( c ) 



1 1 

1 1 1/3 

2 11 

3 3 11/3 

4 5 44/3 

5 12 62 

6 35 814/3 

7 98 3652/3 

8 299 5563 



Notice that one advantage of the triangle-based ex- 
pansion in the triangular lattice, over the square-based 
expansion in the square lattice, is that in the former the 
maximum number of lattice sites of a cluster with N t tri- 
angles is 2N t + 1 , while in the square lattice it is 3N S + 1 
(N s being the number of squares). This means that one 
can fully diagonalize clusters with more triangles than 
squares, which helps both for the bare NLC sums as well 
as for extrapolations. In the next subsections we apply 
the above expansions to Heisenberg and Ising models on 
the triangular lattice. 

A. Heisenberg model 

The triangular-lattice antiferromagnetic Heisenberg 
model is a fascinating quantum spin model, which has 
long-range order at T = [20T | but with spin-spin corre- 
lations that remain short range down to fairly low tem- 



peratures [2l[. In contrast to the square lattice AFHM, 
the AFHM on the triangular lattice shows no evidence 
for a renormalized classical behavior [22, HH, 0] down 
to lowest temperatures that can be reached in HTE. It 
is a frustrated spin model so that QMC methods suffer 
from a sign problem. It has recently been argued [25LI26T] 
that the anomalous finite-temperature behavior is due to 
the excitation of rotons, which leads to high entropy at 
relatively low temperatures. 

The specific heat of the triangular lattice AFHM is a 
quantity that could be of direct experimental interest. 
It has also been calculated from HTE by the recent ap- 
proach of Bcrnu and Misguich (BM) |9fl. These authors 
found that direct Pade approximants [2l[ not only fail 
at surprisingly high temperatures but are also unable to 
correctly locate the maximum of the specific heat and its 
height. 
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FIG. 10: (Color online) Specific heat as a function of temper- 
ature for the Heisenberg model on a triangular lattice. Direct 
sums are compared with BM results in Ref. 0]. 

In Fig. [10] we show the bare sums for the three possi- 
ble NLC expansions discussed before, and compare them 
with the BM results The bond and site expansions, 
up to 12 bonds and 13 sites, respectively, are well con- 
verged only at high temperatures (up to T ~ 2), with 
the site expansion being slightly better. On the other 
hand, the triangle based expansion converges down to a 
lower temperature T ~ 0.6. This temperature is very 
close to the lowest temperature up to which direct Pade 
extrapolations agree with BM results. 

On performing extrapolations over the bare sums 
shown in Fig. [TO] we found that the site expansion is the 
one that enables an improvement of the convergence to 
the lowest temperature. (The bare results for the triangle 
based expansion can hardly be improved by sequence ex- 
trapolation methods.) In Fig. [11] we show results for two 
possible extrapolations of the site based NLC expansion 
compared with BM results @ . 

Notice that we have included two terms for each ex- 
trapolation scheme. To understand what they mean one 
has to realize that up to 13 sites Euler transformation 
allows for 13 terms, from which we have taken the first 
four to be the bare sums, and starting with the fifth we 
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have applied the transformation as explained in Sec. IIHI 
Hence, in Fig. [TT] we are showing the last two terms (in 
Sec. IIVI we showed only the last one). For Wynn's ap- 
proach on the other hand, one should remember that two 
terms are lost after each cycle of improvement. So after, 
five cycles (the case in Fig. ITT]) ten out of the initial 13 
terms in the bare sum are lost, i.e., in Fig. [TTJ we are 
showing the last two of the remaining three. 
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FIG. 11: (Color online) Specific heat as a function of temper- 
ature for the Heisenberg model on a triangular lattice. Ex- 
trapolations are compared with BM results in Ref. [§[. The 
superindex refers to the term in the extrapolation (see text 
for details). 

Figure[TTJshows that while Euler transformation allows 
one to extend the convergence of the site based expan- 
sion to the region where the triangle based expansion was 
convergent, Wynn's extrapolation scheme enables one to 
obtain results at lower temperatures (up to T ~ 0.3). In 
contrast to direct Pade for HTE, Wynn's scheme for NLC 
allows one to reach the maximum of the specific heat as 
predicted by BM @. 

A second quantity of much experimental interest is the 
uniform susceptibility [27j • In Fig. [12] we show NLC (bare 
and extrapolated) results for the uniform susceptibility 
(x) of the AFHM on the triangular lattice. NLC results 
are compared for this quantity with series extrapolations 
of HTE obtained by integrated differential approximants 
[ItJ • (Notice that the BM approach 0] is not suitable for 
calculations of xO 

A comparison between the results for x an d C v helps 
to understand why the flexibility one has for selecting 
different kind of clusters in building the NLC expansion 
can be useful. In contrast to C„, the results of the bare 
sums for the site based and triangle based NLC expan- 
sions for the susceptibility converge well up to about the 
same high temperature (see Fig. [12] vs Fig. [TO]) . This 
might suggest that the triangle based expansion may not 
bring any advantage for this quantity. However, in con- 
trast to C v , series extrapolations extend the region of 
convergence for \ f° r the- triangle based NLC (T Wynn 
in Fig. I12p , and allow one to reach lower temperatures 
than the extrapolations for the site based expansion (S 
Wynn in Fig. [T2"]) . Notice that in the region where NLC 
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FIG. 12: (Color online) Uniform susceptibility as a function 
of temperature for the AFHM on a triangular lattice. NLC 
bare sums are shown for up to 7 and 8 triangles (7 T and 8 
T in the figure), and up to 12 and 13 sites (12 S and 13 S 
in the figure). The corresponding extrapolations for the site 
based (S Wynn) and triangle based (T Wynn) expansions are 
compared with series extrapolation results of HTE, obtained 
by integrated differential approximants [2^]. For the latter, 
only the upper and lower boundaries are shown. 



extrapolations are well converged they are in excellent 
agreement with extrapolations of HTE obtained by in- 
tegrated differential approximants [27| . For x, the inte- 
grated differential approximants of HTE (which do not 
work so well for C v ) appear to have convergence to lower 
temperatures than the ones reached with NLC. 



B. Ising model 

The Ising model [Eq. (fTTI)l on the triangular lattice is 
an exactly soluble model[28. |29| . [30l ] . At zero tempera- 
ture it exhibits power-law decaying spin correlations and 
an extensive entropy S = 0.3231. 
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FIG. 13: (Color online) Entropy as a function of temperature 
for the Ising model on a triangular lattice. NLC bare sums are 
shown for up to 7 and 8 triangles, and up to 12 and 13 sites (12 
S and 13 S in the figure). The corresponding extrapolations 
for the site based (S Wynn) and triangle based (T Wynn) 
expansions are compared with the exact analytical result. 
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For this model, the triangle based NLC expansion con- 
verges down to low enough temperatures that it even 
allows one to study some ground state properties, for ex- 
ample, the entropy shown in Fig. 1131 One can see that 
the site based expansion converges only up to T ~ 1. 
Hence, the triangle based expansion provides one with a 
qualitative improvement over site (and bond) expansion. 
Adding up contributions from clusters up to N t = 1, 2, 
3, 4, 5, 6, 7, 8 triangles leads to estimates of zero tem- 
perature entropy of S = 0.6931, 0.4055, 0.3677, 0.3677, 
0.3499, 0.3521, 0.3417, 0.3440, respectively. In this case, 
the convergence to the thermodynamic limit result is 
power law in 1/N t , compared to an exponential conver- 
gence in the kagome case Q , which is not surprising given 
that the triangular-lattice model is critical [301 ]. 

Wynn's extrapolation of the triangle expansion, up to 
8 triangles, improves towards the thermodynamic limit 
result as shown in Fig. 1131 At low temperatures it gives 
an estimate of the entropy S = 0.3349. The extrapolation 
for the site expansion, on the other hand, only converges 
down to T ~ 0.3. 



TABLE VII: Number of topological clusters and sum of the 
lattice constants for clusters with up to 13 bonds in the 
kagome lattice. The cluster with bonds is the one site graph. 

No. of bonds No. of topological clusters Yl L(c) 






1 


1 


1 


1 


2 


2 


1 


6 


3 


3 


62/3 


4 


4 


77 


5 


8 


304 


6 


17 


3752/3 


7 


36 


5294 


8 


81 


22 845 


9 


194 


299 924/3 


10 


481 


442 507 


11 


1235 


1977 572 


12 


3297 


26 729 935/3 


13 


8944 


40 418174 
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FIG. 14: (Color online) Specific heat as a function of temper- 
ature for the Ising model on a triangular lattice. NLC bare 
sums are shown for up to 7 and 8 triangles, up to 12 and 13 
sites, and compared with the exact analytical result. 

A comparison of the NLC results for the specific heat 
with the exact analytical calculation is shown in Fig. [TJ] 
For this quantity, extrapolations of the site, bond, and 
triangle based expansions do not allow one to improve 
over the direct triangle based sum, so we do not show 
them there. It is remarkable, however, that even though 
the system is critical the results of the triangle based bare 
sums are not far from that exact result. 



VI. KAGOME LATTICE 

In this section we consider the kagome lattice. As be- 
fore, we discuss three different choices of basic clusters 
that allow one to build a consistent NLC expansion. 

The first choice, again, is to use all clusters up to a 
given number of bonds. The number of topological clus- 



ters, and sum of L(c), when grouped by their number of 
bonds is presented in Table. IVIII 

A second choice is to build all clusters with up to a 
given number of sites. When building the Hamiltonian 
for such clusters one needs to place the maximum num- 
ber of bonds possible in them. The number of topologies 
and sum of lattice constants for the site based NLC ex- 
pansion is shown in Table IVIIII We have grouped them 
by number the number of sites. 

TABLE VIII: Number of topological clusters and sum of the 
lattice constants for clusters with up to 15 sites in the kagome 
lattice. 

No. of sites No. of topological clusters ^ L(c) 



1 


1 


1 


2 


1 


2 


3 


2 


14/3 


4 


2 


12 


5 


4 


33 


6 


7 


281/3 


7 


12 


272 


8 


22 


805 


9 


45 


2420 


10 


88 


7358 


11 


183 


22 581 


12 


389 


209 552/3 


13 


842 


217522 


14 


1855 


681 224 


15 


4162 


2 143 905 



Since the kagome lattice consists of corner sharing tri- 
angles, the triangle-based NLC, in this case, involves all 
elementary triangles. This selection of building blocks for 
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NLC reduces dramatically the number of clusters to be 
considered. The number of topologically distinct clusters 
with 1 through 8 triangles, and the sum of their lattice 
constant is shown in Table 11X1 (We have grouped them 
by the number of triangles.) 



TABLE IX: Kagome lattice number of topological clusters 
and sum of the lattice constants for clusters with up to eight 
triangles. The cluster with triangles is the single site. 

No. of triangles No. of topological clusters L(c) 






1 


1 


1 


1 


2/3 


2 


1 


1 


3 


1 


2 


4 


2 


14/3 


5 


2 


12 


6 


5 


94/3 


7 


7 


250/3 


8 


15 


225 



In Ref. [5[ we have already discussed extensively several 
spin models on the kagome lattice. Hence, here we will 
restrict our analysis to the specific heat and the uniform 
susceptibility of the AFHM. 



A. Heisenberg model 

The kagome-lattice AFHM model has been argued to 
have short-ranged spin-spin correlations down to T = 
[31II32I 33]. Its thermod yna mic properties have also been 
of considerable interest [341 ]. In particular, an issue that 
is still under debate (motivated by experiments on He 3 
on graphite) is whether this model exhibits a two-peaked 
structure in the specific heat [35l . l36j . 
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FIG. 15: (Color online) Specific heat as a function of temper- 
ature for the Heisenberg model on a kagome lattice. Direct 
sums, for up to 7 and 8 triangles (7 T and 8 T in the figure) 
and up to 12 and 13 sites (12 S and 13 S in the figure), are 
compared with extrapolations for the triangle expansion and 
with two Pade approximants from Ref. [36(. 



In Ref. [f| we have already studied the specific heat of 
the AFHM. There we compared the direct sums of the 
triangle expansion with Pade approximants from HTE 
(36| . which showed an overall good agreement down to 
T ~ 0.3. The triangle expansion for the specific heat on 
the kagome lattice, in contrast to the triangular lattice 
in the previous section, allows for an acceleration of the 
convergence by means of Wynn's extrapolations. 

In Fig. [T51 we compare the bare sums for the site and 
triangle expansions with results of Wynn's extrapolation 
and Pade approximants [36]. As seen in Fig. [15l the re- 
sults for Wynn's extrapolation appear to converge down 
to T ~ 0.2 and exhibit a clear deviation from the Pade 
results. The deviations are, one could say, in the right di- 
rection since the extrapolation of the specific-heat HTE 
down to T = has a large missing entropy [36| . 
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FIG. 16: (Color online) Uniform susceptibility as a function 
of temperature for the AFHM on a kagome lattice. NLC bare 
sums are shown for up to 7 and 8 triangles (7 T and 8 T in the 
figure), and up to 12 and 13 sites (12 S and 13 S in the figure). 
The corresponding Wynn's extrapolations for triangle based 
expansions are also shown. 

To conclude this section we show in Fig.[l6]NLC results 
for the bare and extrapolated sums of the uniform suscep- 
tibility (x) of the AFHM on the kagome lattice. Similar 
to the extrapolations for the specific heat, Wynn's ex- 
trapolations are well converged down to T ~ 0.2, while 
the direct sums for the triangle expansion converge down 
to T ~ 0.3. On the other hand, the site based expansion, 
up to 13 sites, converges only down to T ~ 1, which is 
the same temperature one can access with HTE with- 
out extrapolations. Overall, for the kagome lattice we 
have found that direct and extrapolated sums of the tri- 
angle based expansion converge better for the thermody- 
namic observables considered (energy, entropy, specific 
heat, and uniform susceptibility) than the site and bond 
based expansions. 



VII. CONCLUSIONS 

We have presented an extensive discussion of the 
numerical linked cluster algorithm introduced in Ref. 
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We have detailed how to construct NLC starting 
from different building blocks on square, triangular, and 
kagome lattices. Several sequence extrapolation tech- 
niques, which we have used to accelerate NLC conver- 
gence, have also been discussed. 

In order to show how NLC works for models with dif- 
ferent ground states and orders, we have studied several 
spin models on square, triangular, and kagome lattices. 
We have shown that NLC is better suited for systems that 
remain short ranged at all temperatures (such as the XY 
model in a staggered field), and for models where corre- 
lations build up slowly so that they become large only at 
very low temperatures. A good example of the latter case 
is the AFHM on the kagome lattice, for which well con- 
verged results could be obtained down to T ~ 0.3 without 
the need of using sequence extrapolation techniques. 

Similar to HTE, NLC also allows for extrapolations be- 
yond the region of convergence provided by clusters up to 
a given system size. It is important to note that within 
NLC the region of convergence is dictated by the largest 
cluster sizes considered, and by the range of correlations 
in the thermodynamic system. Hence, even without ex- 
trapolations, one can, in principle, extend the region of 
convergence by including larger clusters. In this respect 
NLC is fundamentally different from HTE, whose region 
of convergence is dictated by the dominant microscopic 
energy scale, and including larger clusters can only help 
with extrapolations as they do not improve the conver- 
gence of the direct sum. Extrapolations within NLC al- 



low one to go to temperatures lower than accessible by 
means of direct Pade approximants for HTE. Examples 
where NLC is superior in this respect include the specific 
heat in the triangular and kagome lattices. 

Finally, we have also compared NLC results with those 
obtained from exact diagonalization of clusters with pe- 
riodic boundary conditions. We have shown that NLC 
provides accurate results even where ED still suffers from 
very large finite size effects. Even for short ranged mod- 
els such as the XY model in a staggered field, ED can 
fail to predict, for example, the position and height of 
the peak in the specific heat. 

Although it was not implemented here, one way to 
improve NLC convergence at lower temperatures would 
be to use Lanczos type methods to diagonalize larger 
clusters. Larger clusters would become possible if one 
was to focus only on low lying states rather than the full 
diagonalization that we have used in this work. 



Acknowledgments 

This work was supported by the US National Science 
Foundation, Grant Nos. DMR-0240918, DMR-03 12261, 
and PHY-0301052. We are grateful to B. Bernu and G. 
Misguich for providing us with their data from Ref. 
and to R. Yu for providing us with the QMC results pre- 
sented in Sec. [TV] 



[1] J. Jaklic and P. Prelovsek, Adv. Phys. 49, 1 (2000). 

[2] A. W. Sandvik, Phys. Rev. B 59, R14157 (1999). 

[3] For a general introduction see C. Domb and M. S. Green, 
Phase Transitions and Critical Phenomena, (Academic 
Press, New York, 1974), Vol. 3. 

[4] A. J. Guttmann, in Phase Transitions and Critical Phe- 
nomena, edited by C. Domb and J. Lebowitz (Academic 
Press, London, 1989), Vol. 13. 

[5] M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. Lett. 
97, 187202 (2006). 

[6] A similar idea for ground state properties was considered 
earlier by A. C. Irving and C. J. Hamer, Nucl. Phys. B 
230, 361 (1984). 

[7] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. 
T. Vetterling, Numerical Recipes in Fortran (Cambridge 
University Press, Cambridge, England, 1999), Sec. 5.1. 

[8] J. Oitmaa, C. Hamer, and W-H. Zheng, Series Expansion 
Methods for Strongly Interacting Lattice Models (Cam- 
bridge University Press, Cambridge, 2006). 

[9] B. Bernu and G. Misguich, Phys. Rev. B 63, 134409 
(2001). 

[10] R. R. P. Singh, Phys. Rev. B 39, 9760 (1989); J.-K. Kim 

and M. Troyer, Phys. Rev. Lett. 80, 2705 (1998). 
[11] A. W. Sandvik and J. Kurkijarvi, Phys. Rev. B 43, 5950 

(1991); A. W. Sandvik, J. Phys. A 25, 3667 (1992); Phys. 

Rev. B 59, R14157 (1997). 
[12] E. Y. Loh, Jr., J. E. Gubernatis, R. T. Scalettar, S. R. 

White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B 



41, 9301 (1990). 
[13] P. Henelius and A. W. Sandvik, Phys. Rev. B 62, 1102 
(2000). 

[14] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 

(2005) . 

[15] R. Yu (private communication). 

[16] V. G. Rousseau, D. P. Arovas, M. Rigol, F. Hebert, G. G. 
Batrouni, and R. T. Scalettar, Phys. Rev. B 73, 174516 

(2006) . 

[17] A. Priyadarshee, S. Chandrasekharan, J.-W. Lee, and H. 

U. Baranger, Phys. Rev. Lett. 97, 115703 (2006). 
[18] M. Aizenman, E. H. Lieb, R. Seiringer, J. P. Solovej, and 

J. Yngvason, Phys. Rev. A 70, 023612 (2004). 
[19] L. Onsager, Phys. Rev. 65, 117 (1944). 
[20] B. Bernu, C. Lhuillier, and L. Pierre, Phys. Rev. Lett. 

69, 2590 (1992). 
[21] N. Elstner, R. R. P. Singh and A. P. Young, Phys. Rev. 

Lett. 71, 1629 (1993); J. Appl. Phys. 75, 5943 (1994). 
[22] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. 

Rev. Lett. 60,' 1057 (1988). 
[23] P. Azaria, B. Delamotte, and D. Mouhanna, Phys. Rev. 

Lett. 68, 1762 (1992). 
[24] A. V. Chubukov, T. Senthil, and S. Sachdev, Phys. Rev. 

Lett. 72, 2089 (1994). 
[25] W. Zheng, J. O. Fjaerestad, R. R. P. Singh, R. H. McKen- 

zie, and R. Coldea, Phys. Rev. Lett. 96, 057201 (2006); 

Phys. Rev. B 74, 224420 (2006). 
[26] O. A. Starykh, A. V. Chubukov, and A. G. Abanov, Phys. 



14 



Rev. B 74, 180403(R) (2006). 
[27] W. Zheng, R. R. P. Singh, R. H. McKenzie, and R. 

Coldea, Phys. Rev. B 71, 134422 (2005). 
[28] G. H. Wannier, Phys. Rev. 79, 357 (1950). 
[29] R. M. F. Houtappel, Physica (Amsterdan) 16, 425 

(1950). 

[30] J. Stephenson, J. Math. Phys. 5, 1009 (1964); 11, 413 
(1970). 

[31] P. Lecheminant, B. Bernu, C. Lhuillier, L. Pierre, and P. 

Sindzingre, Phys. Rev. B 56, 2521 (1997). 
[32] C. Zeng and V. Elser, Phys. Rev. B 51, 8318 (1995). 
[33] R. R. P. Singh and D. A. Huse, Phys. Rev. Lett. 68, 1766 

(1992). 



[34] P. Sindzingre, G. Misguich, C. Lhuillier, B. Bernu, L. 

Pierre, Ch. Waldtmann, and H.-U. Everts, Phys. Rev. 

Lett. 84, 2953 (2000) 
[35] V. Elser, Phys. Rev. Lett. 62, 2405 (1989). 
[36] N. Elstner and A. P. Young, Phys. Rev. B 50, 6871 

(1994). 

[37] Notice, that from the expression given in Ref. [J] we have 

corrected a sign in a subindex. 
[38] This limitation, however, can be mitigated by the fact 

that one can trivially parallelize the NLC calculation. 

Still, the number of clusters increases exponentially so 

that only few orders can be gained by parallelization. 



